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We study the two-dimensional flow of foams around a circular obstacle within a long channel. 
In experiments, we confine the foam between liquid and glass surfaces. In simulations, we use a 
deterministic software, the Surface Evolver, for bubble details and a stochastic one, the extended 
Potts model, for statistics. We adopt a coherent definition of liquid fraction for all studied systems. 
We vary it in both experiments and simulations, and determine the yield drag of the foam, that is, 
the force exerted on the obstacle by the foam flowing at very low velocity. We find that the yield 
drag is linear over a large range of the ratio of obstacle to bubble size, and is independent of the 
channel width over a large range. Decreasing the liquid fraction, however, strongly increases the 
yield drag; we discuss and interpret this dependence. 



I. INTRODUCTION 

Multiphase materials such as colloids, emulsions, 
polymer or surfactant solutions, wet granular systems 
and suspensions of dcformable objects like red blood 
cells are characterized by a complex mechanical be- 
haviour [l||, due to the interaction of their constitutive 
entities. The concentration is one of the key parameters 
which control the rheology, determining especially the 
transition from liquid-like to solid- like properties [2| ■ 

Amongst these complex fluids, liquid foams provide 
a convenient model experimental system for laboratory 
studies of the interplay between structure, concentra- 
tion and rheology. This is because the bubbles which 
constitute the foam's internal structure can be easily 
visualised and manipulated. The mechanical behaviour 
of foams is very diverse: they appear elastic, plastic 
or viscous depending on the deformation and velocity 
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gradient 0,Q]. 

A liquid foam consists of gas bubbles separated by 
a connected network of liquid boundaries. This liquid 
phase occupies a fraction $ of the volume of the foam. 
The "dry foam" limit, in which <f> tends to zero, corre- 
sponds to polyhedral bubbles separated by thin walls. It 
is associated with a divergence of certain contributions 
to the viscous dissipation [5|. However, the foam's non- 
dissipative properties (such as surface energy shear 
modulus or yield stress 0, Q|) usually tend to a regular, 
finite limit when the liquid fraction (f> tends to zero. 

The total "yield drag" Fy is the minimal force ob- 
served when there exists (or, equivalently, required to 
create) a movement of the foam relative to an obsta- 
cle Q. It is a global, geometry-dependent quantity di- 
rectly measurable in experiments and in practical appli- 
cations of foams, for instance when a foam flows through 
a porous medium [10[, or when one introduces an ob- 
ject into a foam (analogous to sticking one's finger into 
shaving cream). 

In the low-velocity limit (in which viscous dissipation 
is neglected [ill Q) the total yield drag Ft has two 
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bubbles, denoted Fy, and the network of bubble walls 



(i.e.; soap films with surface tension), Fy- Thus 



Fi. = FP 



(1) 



Here we consider the network contribution Fy and show 
how it is affected by the liquid content of the foam. 

We consider a single layer of equal-area bubbles to 
facilitate preparation and analysis of experiments, as 
well as numerical and analytical modelling 13|. Sec- 



tion [IT] presents a 2D flow of a quasi-2D foam (a bub- 
ble monolayer) around a fixed circular obstacle within a 
long channel: this is the historical experiment of Stokes, 
already adapted to foams both in 2D and 3D 

15, Q, 0,01 flows. We compare them with truly 2D 
simulations using two physically equivalent but differ- 
ently optimised software packages (Section IIIip . The 
simulation methods allow easy variation of the geomet- 
rical parameters such as bubble, obstacle and channel 
size and better control of bubble area. In Section llVl we 
discuss the issue of a common, unambiguous definition 
of liquid fraction for all systems, in theory, experiments 
and simulations. Section |V] presents our results: we 
show that the yield drag displays the expected depen- 
dence with the bubble, obstacle and channel size, and 
increases when the liquid fraction <I> decreases. The 
discussion in Section IVII emphasises that taking into 
account the effect of liquid fraction allows all data to 
be plotted on a single master-curve and that, although 
they cover different ranges of <&, the results of both 
simulation and experiment are consistent with a sim- 
ple model. 

II. EXPERIMENTAL METHODS 
A. Foam channel 

Our bulk soap solution is de- ionised water with 1% 
Teepol, a commercial dish-washing liquid. Its surface 
tension, measured with the oscillating bubble method 
(that is, imaging the interface shape), is 7 = 26.1 ± 0.2 
mN m^^, and its kinematic viscosity, measured with a 




FIG. 1: Image of the experiment, with the foam confined 
between liquid and glass and flowing from top to bottom. 
Foam tliickness h = 4.5 mm; bubble area A = \& mm^; 
obstacle diameter do = 3 cm; mean velocity v = 5.6 mm 
s~^; effective liquid fraction <1> — 0.06. 



The experimental set-up [9| confines the foam be- 
tween a liquid reservoir and a glass lid ("liquid-glass" 
set-up [0]). Aim long, Wc = 10 cm wide tank is filled 
with soap solution, leaving below the glass lid a free 
space of thickness h which we can adjust. We will call 
this parameter the "foam thickness" for simplicity. At 
its centre is a circular obstacle of diameter do = 3 (Fig. 
[T|) or 4.8 cm. At the entrance to the channel, nitrogen is 
blown at a computer-controlled flow rate, which varies 
between 5 and 500 ml min^^. A typical value of the 
average velocity is 3 mm s~^, for a 3.5 mm thickness 
and a flow rate of 50 ml min^^. 

The resulting foam consists of a horizontal monolayer 
of bubbles. It exits freely at atmospheric pressure at 
the open end of the channel, P = Patm- In the ab- 
sence of the obstacle, it yields a two-dimensional plug 
flow. With the obstacle present, the flow remains two- 
dimensional (even though the foam itself is not exactly 
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Due to the presence of the obstacle, there is a veloc- 
ity gradient. There arc many bubble rearrangements 
(so called "Tls" or neighbour-swapping events): two 
three-fold vertices contact, merge and re-separate. We 
observe well-separated Tls; that is, between two Tls, 
there is enough time for the foam to relax to an equi- 
librium state. The present flow is slow enough that we 
can extrapolate the results to the low velocity limit , 
where a comparison with quasi-static calculations and 



simulations makes sense 



20|. (Although note that this 



is distinct from the zero-velocity case, that is the ab- 
sence of flow.) The data presented below are all in this 
limit, extracted from the experiments as described in 
Appendix lA II 

The bubble walls meet the solid boundaries of the 
foam (glass plate, lateral channel walls, obstacle itself) 
at a 90° angle 



2l| . The surface density of bubbles is 
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23| is thus negli- 



1/A, where A is the average area per bubble (including 
its walls). The foam is monodisperse: the area variation 
at the channel entrance is less than 5%. The average 
area is fixed at a value ranging from 0.121 cm^ to 0.393 
cm^; most experiments have A = 0.160 cm^. Despite 
the low velocity, and hence the long transit time, we 
detect neither bubble coalescence nor coarsening. The 
effect of foam ageing on rheology 
gible. 



B. Force measurements 

1. Total yield drag 

The obstacle floats just below the top glass surface 
and is free to move, without solid friction. However, it is 
linked to a fixed base through a calibrated elastic fibre. 
We track the obstacle displacement from its position at 
rest using a CCD camera which images the foam flow 
from above. We thus measure the force exerted by the 
flowing foam on the obstacle (precision better than 0.1 
mN) [g. 

We check that the lift (spanwise component of the 



as expected by symmetry (data not shown). After a 
transient, the total drag (strcamwise component of 
the resultant force) fluctuates around a steady value: 
we record the average and standard deviation of these 
steady flow data. The extrapolation to the low velocity 
limit (or zero-velocity intercept) of the force-velocity 
curve deflnes the yield drag Ft- It is independent of 
the bulk solution viscosity [24 , and increases with the 
obstacle to bubble size ratio [9|. 

In this paper, we reanalyse the data already pub- 
lished in [9[ at various bubble areas, and we present 
new data for another control parameter: the foam thick- 
ness. These data are presented in Appendix lA II As 
explained in Sec. IIVB| the foam thickness provides a 
means by which to vary the liquid fraction in experi- 
ments. 



2. Network contribution to the yield drag 

We measure Fy as follows. Each bubble wall in con- 
tact with the obstacle pulls it with a force equal to its 
line tension A (the energy per unit length, which is of or- 
der 2^h, see Appendix I A 2p . The elastic contribution of 
the wall network to the drag is then the vectorial sum of 
all these individual forces, which all have the same mod- 
ulus A. As mentioned above, in a quasi-static flow each 
wall touches the obstacle at 90° angle. Thus it suffices 
to find the contact points between bubble walls and the 
obstacle and sum all outward normal vectors to the ob- 
stacle vectorially at these contact points (which is easy 
to determine for a circular obstacle) . If the downstream 
geometry of the foam was the same as that upstream, 
the drag would be zero. Since bubbles are squashed 
upstream and stretched downstream, the asymmetry 
means that there are more bubble walls pulling the ob- 
stacle downstream, and we measure a downstream elas- 
tic contribution to the drag, Fy/X. 

The actual value of the line tension A is unimportant 
in what follows, where only measurements of Fy/X are 
compared. However, as presented in detail in Appendix 
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of magnitudes of the independent measurements of Fy 



and F^. 



III. SIMULATIONS 

A. Deterministic simulations (the Surface 
Evolver) 



The Surface Evolver 



25 



26[ offers the possibility to 



reach a true quasi-static limit, that is a succession of 
exact equilibrium states, through a deterministic min- 
imisation of the foam's energy. It yields precise details 
of the foam structure. 



(a) 




1. Preparation of the foam 

We use a mode in which all bubble walls are rep- 
resented as circular arcs. The Surface Evolver lets 
these circular arcs evolve in order to minimise the to- 
tal perimeter (equivalent to the energy, up to the pref- 
actor A). It enforces the constraint that bubble areas 
A remain fixed and determines the corresponding La- 
grange multipliers, namely each bubble's pressure P. 
Since we can freely choose the units, we call them "cm" 
and use bubble size A — 0.16 or 0.353 cm^. channel 
width Wc ~ W cm, and obstacle diameters do = 1.5, 3 
and 4.8 cm, to reproduce actual experiments. 

The lateral sides of the channel are rigid and do not 
interact with the foam, ensuring free-slip boundary con- 
ditions for the flow, resulting in a 90° angle where a 
bubble wall meets the side. We adopt a periodic bound- 
ary condition in the direction of motion: bubbles that 
exit at the end of the channel are fed back into the en- 
trance of the channel. We stop the simulation when 
each bubble has passed the obstacle no more than once. 

We begin with a rectangular lattice of 30 x 25 
monodisperse bubbles of area slightly larger than the 
required area A. We randomly perturb this lattice so 
that all the unstable four-fold vertices dissociate into 
pairs of three-fold vertices and the whole foam struc- 



(b) 




FIG. 2: Images of simulated foam flow. The x-axis is par- 
allel to the flow along the channel, with periodic bound- 
ary conditions (exiting bubbles re-enter); axis y is spanwise, 
with free-slip rigid boundary conditions on either side of 
the channel, (a) Surface Evolver. The image shows the 
whole simulation domain of 750 bubbles; the shaded bub- 
bles started in a horizontal line. Here do = 4.8 cm, A = 0.16 
cm'^, Wc — 10 cm, Lc = 0.05 cm and therefore $ = 0.0037. 
(b) Potts model. The image shows the simulated channel's 
full width (except for a few pixels) of 256 pix, and half its 
length. Here do = 74 pix, A = 100 pix^ and $ = 0.005. 
Bubbles coloured in white are without topological defect: 
6-sided bulk bubbles, or 5-sided bubbles touching a lateral 
wall or the obstacle [g]. Bubbles with fewer neighbours are 
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bubble to be the circular obstacle, and slowly increase 
its area to the required value (and correspondingly re- 
duce the bubble areas to A) and constrain its edges to 
lie on a circle. The centre of the circular obstacle is 
then moved to the centre of the channel and the struc- 
ture again relaxed to equilibrium. 

2. Simulation of the flow 

With the obstacle in the desired location and the 
foam close to equilibrium, we start the quasi-static iter- 
ation procedure. This requires that we move the foam 
past the obstacle, in a direction which we denote by x] 
the difficulty is in doing this with the periodic boundary 
conditions without fixing any vertices or bubble shapes. 
Our method is to choose a continuous line of consecu- 
tive bubble walls from one side of the channel to the 
other. Joining this to a line at a; = with lines along 
the channel walls defines a plane region with a certain 
area, which we constrain. At each iteration we choose 
a convenient line of consecutive walls, and increment 
the target area of the region formed by a small amount 
dA (equal to 0.05 cm^ in all simulations), resulting in 
a slight movement of a line of films without modifica- 
tions to the bubble areas. The total perimeter of the 
structure is then reduced until it converges to a con- 
stant value (Fig. [5^), so that measurements can be 
performed. 

We have double-precision values for the network ge- 
ometry. We measure the network contribution Fy to 
the yield drag as in experiments (section III B 2|1 . It is 
the sum of the unit vectors of the bubble wall with one 
end attached to the obstacle, expressed in units of the 
line tension (hence as a dimensionless number). Here 
too, we check that the lift is consistently zero within 
fluctuations (data not shown). 

With the area increment dA = 0.05 cm^, the tran- 
sient lasts for about 600 iteration steps (Fig. [3^). This 
is comparable to, but still smaller than, the total simula- 
tion time that is reasonably accessible. After this tran- 



(a) 




400 600 800 1000 
Number of iterations 



(b) 



0) 



o 

0) 

2; 



t,'L,l!"'l 



Hi 



100000 200000 300000 400000 500000 600000 
Monte Carlo Steps 

FIG. 3: The network force Fy (expressed in units of the 
line tension A) measured in simulations versus time, (a) 
Surface Evolver data plotted every iteration step; do = 4.8 
cm, A = 0.16 cm^ and Lc = 0.05 cm, $ = 0.004. The 
plateau value is Fy ~ 9.6±1.4. (b) Potts model data plotted 
every 1500 Monte Carlo Steps; do = 74 pix, A = 100 pix^, 
$ = 0.005. The plateau value is F{^ = 5.4 ± 1.1. 



fluctuations, due to the rearrangements of the bubbles, 
recall the stress drops observed in Couette experiments 
for disordered foams 0, [28 1. We record the average 
and standard deviation of these plateau (steady-flow) 
data for a total of 1500-600 = 900 iterations. To vali- 
date the choice of our simulation size, we checked once 
that the drag forces are the same with more bubbles 
in the direction of flow (1250 bubbles instead of 750), 
although the transient is longer. 

Each simulation takes about 35 hours on a Pentium 
IV 3.20 GHz processor: typically a several hour build- 
up to the initial structure (inflating the obstacle), plus 
one iteration per minute (depending on the number of 
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B. Stochastic simulations (Potts Model) 



To simulate a larger number of bubbles, the Potts 



model adapted for foam rheology [29| also minimises the 
same energy, but stochastically (Monte-Carlo), which 
increases the simulation speed. It thus provides more 
statistics on Fy and allows quicker variation of the ge- 
ometrical parameters. 

1. Principle of the Potts Model 

The Potts model is derived from a large-Q Potts 
model run at zero temperature, a model widely used 
to model grains in crystals 3^. It has been also ap- 
plied to different domains of foam physics, including 
rheology, by enforcing the conservation of bubble size 
and applying an external force [3| • 

We consider a 2D square lattice. Each site i has an 
integer index (7^. The k*'^ bubble is defined as the do- 
main consisting of all sites with the same index value 
Gi = k. Thus bubbles tile the plane without gaps or 
overlaps. The evolution is driven by the minimisation 
of a total energy 7i (strictly speaking, it is a Hamil- 
tonian), which has the same three physical ingredients 
as in the Surface Evolver: interfacial energy, area con- 
straints, external forcing of the flow. Since the calcula- 
tions are performed on a lattice, we have 

n = \ [i-'5(a„a,)] 

i,j neighbours 



+x E 



A 



_ (2) 

bubbles k sites i 

The first term represents the contribution of the en- 
ergy of the interfaces between the bubbles. Minimising 
this term leads to perimeter minimisation. Here 5 is the 
Kronecker symbol: 1 — (5 is equal to 1 if the neighbour- 
ing sites i, j belong to different bubbles (ci ^ aj); else 
it equals zero. We choose to evaluate this term with 
the fourth nearest neighbour interactions to obtain an 
isotropic line tension insensitive to the details of the 
lattice [sil. 

The second term keeps each bubble Ak (the num- 



target value Aj^. Here x is the compressibility, which 
we choose to be high enough to keep bubble areas con- 
stant to within a few pixels. The balance between this 
term and the preceding one simulates a foam relaxing 
towards mechanical equilibrium. 

The third term is a bias term that describes an en- 
ergy gradient, hence a homogeneous external force field. 
Here b is the bias intensity and x the site's coordinate 
along the flow. Without obstacle, the resulting velocity 
profile would be a plug flow. 

We use a Metropolis algorithm to evolve the foam: 
we randomly select a site at a bubble boundary, change 
its index to the value of a neighbour if and only if this 
decreases the total energy (cq. [2]). Several indepen- 
dent changes are tried successively; a Monte Carlo Step 
(MCS) is defined conventionally as a number of tries 
equal to the total number of lattice sites. 

2. Simulation of the flow 



As for the Surface Evolver (section rill A we choose 
a periodic boundary condition in the direction of flow 
and free-slip rigid boundary conditions on the channel 
sides. To ensure that it docs not affect the steady-state 
measurements presented below, the total channel length 
is 4u'c, out of which only 2'Wc arc used for measurements 
and are shown on Fig. ([lb)- 

To match the experiments, we choose 16 < do < 148 
pix, 64 < A < 400 pix^ and 64 < < 512 pix. Ini- 
tially, we insert a rigid round obstacle in the centre 
of the channel, and let a perfectly ordered foam (hon- 
eycomb pattern) flow in. We then switch off the bias 
term by setting 6 = 0, and relax the foam to ensure that 
the bubbles recover their (near) equilibrium state. The 
foam has reached the stationary state at the end of this 
preparation. 

We then switch the bias on again, and use the small- 
est bias b for which the foam flows, which is constant 
and independent of parameters such as bubble diam- 
eter. We perform measurements at intervals of 1500 
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We measure the network contribution to the drag us- 
ing the same method as in the experiments and Sur- 
face Evolvcr simulations. It fluctuates around a steady 
value: we record the average and standard deviation of 
these plateau (steady-flow) data (Fig. ^jp) . We run each 
simulation for a total of 600,000 MCS, during which a 
bubble passes completely through the channel but no 
bubble passes the obstacle twice. One simulation takes 
about 12 hours on a Pentium IV 2.8 GHz processor. 

IV. LIQUID FRACTION 



In ideal 2D foams (Sec. IIV Ap . given A, the liquid 
fraction scales as the square of the vertex radius. This 
radius in turn relates to a cut-off length Lc, that is the 
length at which a bubble edge becomes unstable and the 
surrounding bubbles undergo a Tl. This length Lc is 
always defined, and is relevant to the mechanical prop- 
erties investigated here. In experiments, Lc depends on 
the actual (3D) shape of bubbles; in simulations, Lc is 
an input parameter. 

It suggests a coherent definition of the effective liq- 
uid fraction, presented below, consistent within and be- 
tween experiments, simulations and theory. Note that 
we consider here a foam dry enough to have a non-zero 
shear modulus 32], that is, below the critical liquid 



fraction Q (see eq. [8] for the honeycomb value). 

A. Ideal 2D foams 

In an ideal 2D foam, a Plateau border is a triangle 
with concave edges of radius R which match tangen- 
tially three straight lines meeting at 120° (Fig. |4|). The 
area Aps of a Plateau border is 



The liquid fraction $ is defined as: 
A, = A^, 



(3) 



(4) 




(b) 




FIG. 4: Picture of two adjacent three-fold vertices with 
Plateau borders, (a) We apply the decoration theorem 
[s^ l to model a wet foam. The liquid is present only at 
the vertices, and (assuming here straight or nearly-straight 
walls) the uniformity of pressure P inside bubbles forces 
each gas/liquid interface to have the same radius of curva- 
ture, R. (b) Critical position of the vertices just before the 
"Tl" neighbour-swapping event. It defines the cut-ofF wall 
length Lc- 

tween 3 bubbles, we have Ai = nABp/i- For a honey- 
comb array of bubbles, n = 6 and: 



(2V3^ 



i?2 

A ' 



(5) 



where Ai is the area occupied by the liquid. Since each 



Here we consider foams where bubbles have the same 
area (monodisperse foams), but not necessarily the 
same number of sides n (topological disorder). Since 
on average over the whole foam n « 6 and since 
Plateau borders have almost the same size and radius 
of curvature, eq. ^ still holds here approximately. 

A Tl is triggered when the distance between these 
vertices becomes smaller than a cut-off wall length Lc, 
which increases with R, and thus with $. To make this 
observation more quantitative, one possible convention 



8 



(Fig. ED: 



the systematic error on Lmax- The total uncertainty on 



^/3 2 ' 



so that, together with eq. 



$=3 /V3-^) 4^ 0.242 ^. 
2 V 2/ A A 



(6) 



(7) 



Given A, the physical information conveyed by R, $ 
or Lc is the same. For comparison between different 
experiments or simulations, we use $ because it is di- 
mensionless. 

For an ideal honeycomb without shear, all vertices 
merge at the same liquid fraction: the hexagons become 
circular when Lc equals the side-length of the hexagons; 
the bubbles are circular, with a radius equal to R. This 
critical liquid fraction is Q: 



2V3 



0.0931. 



(8) 



B. Experiments 

In the experiments, the actual (3D) shape of bubbles 
is determined by the foam thickness. Fig. ([T]) shows 
a foam thickness of 4.5 mm; beyond this thickness, the 
bubbles undergo a three-dimensional instability and the 
foam is no longer a monolayer 3^. At the other ex- 
treme, at 2 mm and below, the bubbles are circular and 
separated: both the foam's 2D shear modulus and the 
yield drag vanish. When h increases, Lc decreases, thus 
$ decreases too (Fig. [9] in Appendix . 

We show in Appendix |B] that there is a correspon- 
dence between Lc and the length Lmax of the edge of a 
bubble just attached to the obstacle (Fig. [TTb in Ap- 
pendix |B|. Since Lmax is much bigger than L^ it can 
be estimated with much more precision in experiment 
(the uncertainty in both Lc and L^ax is one pixel). 

We measure on the skeletonized image the length 
imax several times preceding its disappearance during a 
Tl process, and keep the average as Lmax and the stan- 
dard deviation as SLf^^. The skeletonization itself in- 
duces a systematic error in determining the actual posi- 



Lniax is therefore (5Lmax 



{SL^Ilr + ('5imax 



skol ^2 



To deduce $ from the measurements of Lmax, we com- 
bine Eqs. ^ and (jB3P to get the following expression 
as a function of L'^^^^/A only: 



3 2^/3- 



2 2 + x/3 \Lr, 



Its uncertainty is: 

^ =2 ^ 
$ L^ 



A/L, 



4V3A 



imax/4V3A 



(9) 



Va/l„ 



imax/4V3A 



C. Simulations 

The Surface Evolver requires that we specify explic- 
itly the cut-off wall length Lc at which two three-fold 
vertices are allowed to contact, merge and re-separate. 
Since Lc is an input parameter, it is determined with- 
out uncertainty. This defines explicitly an effective liq- 
uid fraction (eq. [7]), at least for small values of $. We 
choose Lc to be of the order of 0.1 cm or slightly smaller, 
reaching $ = 0.0015, 0.0037, 0.0061 and 0.015. At very 
small values of $ < 6 10~^, films behind the obstacle 
would get very stretched and lead to numerical prob- 
lems. Attempting larger values of $ > 0.015 would 
lead to poor convergence in the Surface Evolver and 
would require that we simulate the actual geometry of 
the liquid in the vertices (including 4-fold vertices). 

In the Potts model, the cut-off distance Lc at which 
two vertices merge is either 1 or 2 pixels. We use this 
range to define the uncertainty on the value of <&. Since 
we will plot the results in log scale, we choose Lc ~ \/2 
and $ w 0.242 x 2/A w 0.5^-^ to he in the middle 
of this interval. Thus the area A of bubbles (that is, 
the number of pixels per bubble) defines an effective 
liquid fraction, at least for small values of $. The sim- 
ulated range 64 pix^ < A < 400 pix^ corresponds to 
0.00125 < <f> < 0.0075, large enough to describe realis- 
tically the shape of bubbles, and small enough to keep 
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RESULTS 



The experiments and both simulations present qual- 
itatively similar images (Figs. [1] ^ and consistent re- 
sults for the yield drag force, always directed down- 
stream. 

There are a priori four lengths in this problem: the 
channel width Wc, the obstacle diameter do, the bubble 
size ^/A; and the cut-off length Lc- As far as we can 
tell, it is safe to assume that the channel length (if long 
enough) is irrelevant here. These four lengths can be 
reduced to three dimensionless parameters. We present 
the results using: do/wc which characterises the flow 
geometry; do/^/A which describes the foam-obstacle in- 
teraction; and L^/A which characterises the threshold 
for Tl rearrangements, and corresponds to the liquid 
fraction $. 

A. Effect of obstacle to channel size ratio 

Potts model simulations indicate that the network 
yield drag is independent of the ratio of obstacle size to 
channel width, d^/wc (Fig. [5^). This ceases to be valid 
at small do, when the obstacle is comparable in size to 
a bubble, and at large do, when the distance between 
the obstacle and the channel side is small [2^. 

The lack of dependence on do / Wc that we find char- 
acterises the yielding behaviour of the foam: it means 
that only a small region near the obstacle is affected by 
the flow Q|. Nonetheless, the zone where the obstacle 

34 1 than in 3D Q, 



influences the flow is larger in 2D 



as elastic or hydrodynamic interactions would suggest. 

B. Effect of obstacle to bubble size ratio 

Potts model simulations indicate that the network 
yield drag increases linearly with the obstacle size do 
(Fig. Eb) at fixed bubble area. This is consistent with 
the force increasing as do/v^, also suggested by the 
available Surface Evolver data, as well as by experimen- 
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FIG. 5; Network contribution to the yield drag Fy (ex- 
pressed in unit of A), measured in Potts model simulations 
with A = 100 pix^ ($ = 0.005). (a) versus Wc, for 
do = 74 pix; the solid line is the average (value 5.43). (b) 
Fy versus do, for Wc — 256 pix; the solid line is a linear fit 
with zero intercept, Fy =0.77 do/\/3. 

of the obstacle's spanwise dimension ("leading edge") 
Note that most elastic properties of a foam scale 
like 1/VA In fact, when A increases, the density 
of bubbles and of bubble walls decreases, and so does 
a foam's elastic modulus (it would eventually vanish if 
there were only one large bubble left). 



C. Effect of liquid fraction 



We need to separate the effects of foam geometry. 
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Effective liquid fraction 
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LllA 

FIG. 6: The network contribution to tfie yield drag Fy, 
rescaled by Ac?o/vC4, is plotted versus the dimensionless 
quantity L'^/A (bottom scale), that is versus the effective 
liquid fraction $ (top scale), in the range 10^"^ < $ < 10^^. 
All control parameters (do, Wc, A and <1?) are varied. Ver- 
tical bars indicate the standard deviation of the force fluc- 
tuations in time around the plateau value. Horizontal bars 
indicate the uncertainty discussed in Sec. IIVI Data are 
from experiments (■), Surface Evolver (o) and Potts model 
(a); the solid line denotes the analytical model, from eq. 
(|B10[) . without adjustable parameters. Note that the hori- 
zontal scales are logarithmic and shifted with respect to each 
other. 

rescale the network contribution to the yield drag Fp 
by cIq/^/A, and plot all our data as a function of $. All 
the data, from both experiments and simulations, are 
well rescaled in the range lO"'"* < <f> < 10^^ (Fig. O. 
This is the main result of the present paper. 



VI. DISCUSSION 
A. Model 

The effect of the liquid fraction on the network drag 
can be understood as follows. A bubble of area A de- 
taches from the obstacle when its width is of order Lc, 
and thus its length is of order A/Lc- When <I> decreases, 
Lc decreases too. Bubbles stretch more downstream. 



ber of bubble walls pulling the obstacle downstream in- 
creases; simultaneously, the number of walls upstream 
decreases. This larger up/downstream asymmetry re- 
sults in an increase in the resulting drag Fy ■ The contri- 
bution from the network (or bubble walls) increases as 
their number per unit length along the obstacle bound- 
ary, namely , and thus scales like 

However, the length of the region on which stretched 
bubbles act decreases, and the divergence in is in 

fact softened by a geometrical factor. As shown in Ap- 
pendix [B] we can estimate this factor by integrating the 
bubble wall contribution around the obstacle. When $ 
increases, Fy decreases; it vanishes for $ = 0.086. This 
is close to the rigidity loss value (eq. [8|). Eq. (|B10|) 
is plotted in Fig. ([S]), without adjustable parameters. 
It shows qualitative agreement with the data over two 
decades of liquid fraction, suggesting that it captures 
the essence of the physics. 



B. Influence of the control parameters 

In the limit of low $, the development of the above 
argument indicates that, provided that the obstacle di- 
ameter and the obstacle-wall distances are larger than 
the bubble diameter, Fy increases according to: 



" $1/4 



0.516 Ado 



(10) 



In simulations, if we multiply the bubble and obstacle 
diameters, expressed in units of the cut-off length, by 
the same prefactor, the network drag changes (data not 
shown), due to the change in $. 

Conversely, increasing only the bubble area A at fixed 
Lc simultaneously decreases both do/\/A and <f>. This 
has two opposing effects, the former decreasing Fy, the 
latter increasing it, resulting in an almost constant Fy 
(Fig. [7]). (In fact there is a weak dependence on area, 
varying as A~^/'^.) This shows that the relevant way 
to vary the liquid fraction in simulations is to modify 
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FIG. 7: Opposite effects in simulations of Fy (here with 
Potts model). When A/L^ increases, both d(,/\/~A and $ ~ 
decrease; so that Fy barely varies (eq. I10|) . Obstacle 
diameter do equal to 16 (a), 32 (♦), 74 (□) and 128 (•). 




Surface Evolvcr simulations allow us to probe the 
range 10^^ < $ < 10^"^. They indicate that the force 
saturates below $ ~ 10^"^, in agreement with our pre- 
liminary experiments of a foam confined between glass 
plates (data not shown). 

Direct observation of simulation images of very dry 
foams (Fig. [51) confirms that the up/downstream asym- 
metry in the number of bubbles touching the obstacle 
is around 10, roughly independent of liquid fraction. 

The model seems to correctly describe the squashing 
and stretching of bubble shapes. However, the inter- 
polation between both extreme values assumes a phe- 
nomenological expression (eq. IB4[) . It seems approxi- 
mately valid only for 10^'^ < $ < 10~^ (see Appendix 
IB]) . It applies to other obstacle shapes, such as an el- 
lipse [35 1 . 




FIG. 8: Surface Evolver simulation of a very dry foam: 
zoom around the obstacle, (a) $ = 3.7 10"'^; (b) $ — 
6 10"®. Here do = 4.8 cm, A = 0.16 cm^, Wc = 10 cm. Flow 
from top to bottom. 



D. Yield drag versus yield stress 



Princen and Kiss [36| have shown that in three di- 
mensions a foam's yield stress scales as: ay — 7(1 — 
'^zdY^^Y {(^-ijj) / R^2^ where i?32 is the surface-volume 
mean radius (Sauter radius), and y(<&3_D) a decreas- 
ing function of <I>3£) which is approximately i^(<&3_D) — 
-0.080- 0.1141n$3£, 136|. 



At this stage, it is worth discussing the fundamental 
differences between yield stress and yield drag. 

The yield stress or yield strain is an intrinsic property 
of the foam. On the other hand, the yield drag depends 
on the geometry of the flow: the foam does not yield 
everywhere around the obstacle, and especially not at 
angles \6\ w 7r/4 from the downstream direction, as ap- 
pears both in experiments (Fig. [1]) and in simulations 
(Figs. [2] and [8]). This spatial dependence implies that 
the relation between yield stress and yield drag is non- 
trivial and ^-dependent. In particular, in a dryer foam 
(Fig. [5|) , the region where bubbles reach their maximal 
deformation is narrower 34 1. 
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only shear, but also elongation, especially near the front 
and back of the obstacle. When $ decreases, the bubble 
elongation can become arbitrary large and dominates 
the contribution to the yield drag. 

VII. CONCLUSION 

To summarise, we investigate the two-dimensional 
flow of a foam around a circular obstacle, within a 
long channel. Our deterministic (Surface Evolver) and 
stochastic (Potts model) simulations, as well as our 
model and experiments, complement and validate each 
other. 

The yield drag is defined as the low-velocity limit of 
the interaction force between an obstacle and a flow- 
ing foam. The network contribution scales as the ratio 
of obstacle to bubble diameter, as long as this ratio 
is larger than unity, and is almost independent of the 
channel width. It increases (because more and more 
stretched bubbles accumulate behind the obstacle) as 
a power law when the liquid fraction contained in the 
foam decreases to 10"'^, then saturates. 

Having found a relevant definition of the liquid frac- 
tion, which is appropriate for experiments, simulations 
and theory, the dependence of yield drag with liquid 
fraction is well characterized. It is very different from 
that of local intrinsic properties such as the yield stress 
or shear modulus. This observation suggests that it 
will be difficult to deduce one quantity from the other. 
This should be kept in mind in future simulations, and 
has to be taken into account when modelling the foam 
behaviour. 

Note that this definition of liquid fraction can be ex- 
tended to other 2D flows in experiments (quasi-2D foam 
set-ups [3|) or simulations. Extension to 3D [l^ 
should also be simple, especially since the main effect of 
quasi-2D set-ups - external friction on the glass plate 



fraction of water |12l|). In simiflations too, our deflnition 
immediately extends to 3D for both Surface Evolver, 
which uses as input parameter the cut-off area for a 
face which undergoes a Tl; and Potts model, where the 
voxel (3D pixel) size plays the same role. 
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- does not seem dominant here. Our present 
effective liquid fraction based on rheological properties 
(Tls) facilitates the comparison between 2D and 3D 



APPENDIX A: FORCE MEASUREMENTS IN 
EXPERIMENTS 

1. Variation with foam thickness 

We present here new data concerning the drag ex- 
erted by a flowing foam of bubble area A ~ 16.0 mm^ 
on a circular obstacle of diameter cIq = 3 cm. We mea- 
sured the drag, as explained in full detail in ^, versus 
the foam velocity V for six different foam thicknesses 
(Fig. [9K). As usual, the drag increases with increasing 
foam velocity. 

More importantly for this paper, at given velocity, 
and especially at the limit of vanishing velocity, the 
drag increases with increasing foam thickness. This is 
due to (i) the decrease of liquid fraction with increasing 
foam thickness, as shown by the snapshots of the two 
extreme foam thicknesses in Fig. (H^); (ii) the increase 
in the height of the films with increasing foam thickness. 
We fit the data by the formula F = F^+AV to get the 
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Velocity (cm/s) 
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FIG. 9: (a) Drag versus velocity for various distances h 
between the top plate and the bottom solution: 2.0 (•), 2.5 
(□), 3.0 (o), 3.5 (A), 4.0 (T) and 4.5 mm (x). Here A = 16.0 
mm^ and do = 3 cm. The curves are the best fits by the 
formula F — Fy + A x . Snapshots of the foam of the 
smallest (bottom) and highest (top) h are also displayed, 
(b) Photograph of a single soap film pulling on the left side 
of the circular obstacle. 

thicknesses. 

Notably, for the smallest foam thickness {h — 2.0 
mm) the foam is almost decompacted and the drag 
tends to vanish at low velocity. More precisely, the 
fit gives an unphysical negative _value. This suggests 
that the rigidity loss transition 



32| occurs for a foam 



thickness between 2.0 and 2.5 mm. 

2. Comparison of network and total yield drags 

We measure the line tension directly as the force ex- 
erted on the obstacle by a single soap film, as shown 
in Fig. ([9)3), for two thicknesses. Its value is 0.44 niN 




FIG. 10: Shape of the interfaces between bubbles, for bub- 
bles of area 16.0 mm^ and volume V — 16.0 x 3.5 mm^, cal- 
culated with the Surface Evolver. Vertical films are shown in 
light grey and the liquid surface in dark grey. To reproduce 
the experiment, we enforce the hexagonal symmetry and in- 
clude the buoyancy. For simplicity, the junction between 
lateral faces and the top plate is assumed to be orthogonal. 



suggests that \/h « 110 mN/m. 

The actual gas-liquid interfaces have a 3D curva- 
ture to match tangentially the water surface and the 
glass plate (Fig. fTU]) . This explains why the the mea- 
sured value is between a lower bound, X/h = 27 « 52 
mN/m expected for a vertical soap film (that is, two 
flat gas/liquid interfaces), and an upper bound, X/h = 
(2 + 7r)7 w 134 mN/m, expected for two films with cir- 



cular cross section (see 39[ for details) 



Using the measured value of A, we determine the ab- 
solute value of the network yield drag Fy in experi- 
ments. We check (data not shown) that Fy is consis- 
tently of the same order of magnitude, but lower than, 
the value of Fy measured directly. The remaining part 
is attributed to the pressure contribution Fy , to be de- 
scribed in a further paper |34l |. The spatial variation 



of bubble height h due to pressure differences is always 
less than 10%, giving an upper limit to the spatial vari- 
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APPENDIX B: VARIATION OF WITH THE 
CUT-OFF LENGTH 

We consider here only the bubbles touching the obsta- 
cle, and the contribution of their walls to the yield drag. 
We assume that (i) the foam is truly 2D and all bub- 
bles have the same area (ideal 2D monodisperse foam) ; 
(ii) pressure P is the same for each bubble, i.e. bubble 
walls are straight and all Plateau borders have the same 
radius of curvature R; (iii) the obstacle is much larger 
than the bubbles {VA <C do) so that we can neglect 
its curvature at the bubble scale. This latter approxi- 
mation could in principle affect the bubbles upstream, 
which share a long edge with the obstacle. However, 
it should not greatly affect the bubbles downstream, 
which are the main contributors to the drag. 



1. Geometry 



(a) 



R 



A 



(b) 



(c) 




2R 



FIG. 11: Model configuration of the bubbles in contact with 
the obstacle, (a) Equilibrium configuration: the dashed line 
represents a polygonal bubble at $ = 0, the dotted line is a 
circular bubble at <I> = <I>c, and the plain line represents the 
intermediate case (0 < $ < "l"c), with straight edges and 
curved triangular vertices, (b) Configuration at the limit 
of detachment: two neighbouring vertices on the boundary 
of the obstacle come into contact, (c) Configuration at the 
point of attachment of a new bubble. There is one vertex 
between bubble 2, bubble 3 and the obstacle, and a second 
vertex between bubbles 1, 2 and 3. When these two vertices 
come in contact, bubble 1 attaches to the obstacle. 



To model a wet foam, we apply the decoration theo- 
rem 32| : the liquid is present only at the vertices which 



decorate an ideally dry foam. For a bubble touching 
the obstacle, we denote by L the distance between two 
neighbouring vertices in contact with the obstacle (Fig. 

m). 

When the foam flows, bubbles attach to the obsta- 
cle upstream, and detach from it downstream. Visual 
observation of both experiments (Fig. [1]) and simula- 
tions (Fig. [2|) indicate that bubbles are flattened along 
the obstacle at the leading side of the obstacle, and 
that they progressively stretch streamwise at the trail- 
ing side. 

L reaches its minimum value downstream, where bub- 
bles detach. There, two neighbouring (decorated) ver- 
tices come in contact, and L equals the cut-off length 
2R (Fig. [ITJd). 

On the other hand, for a new bubble to attach to the 
obstacle upstream, two bubbles must detach through 
the configuration of Fig. (jllb ). In this case, a vertex 
between three bubbles merges with one between two 



and now equals (1 + l/\/3)i?. This geometrically deter- 
mines that the maximum bubble width Lmax obeys: 



Inverting eq. (|Bip yields L 



(Bl) 



imax(A^) = 2^/(^/3 +l)2i?2+Ax/3 

-2{V3 + l)R. (B2) 
At low liquid fraction, imax tends to a finite value. 



namely \/4A\/3; there is no singularity at vanishing R. 
Conversely, at high liquid fraction, L,„ax varies greatly 
with R, so it is preferable to rewrite eq. (|B2p and de- 
termine R from the measurement of imax- 



1 



1 



V3 



A 



in 



4:V3 J 



(B3) 



2. Continuous assumption 

We assume that the shape of the bubbles varies 
smoothly from the configuration of Fig. (jllb) upstream 
to that of Fig. pTb ) downstream: 2R < L < Lmax- 
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switch from a discrete to a continuous description of the 
bubbles. We thus consider L as a continuous function 
of the ortho-radial angle 9 along the obstacle boundary: 
9 = downstream, tt (and — tt) upstream. Equivalently, 
L^^ is the linear density of vertices along the obstacle 
boundary. 

Then L(0) = 2R, L{±t:) = L^ax- To interpolate 
between these values, we assume the following phe- 
nomenological dependence, reflecting that all bubbles 
in the range \9\ > 7r/2 appear squashed against the ob- 
stacle: 



\9\<Tr/2: L{9) = [ R + ^Ln,^.^ 



R - I cos 



\9\>tt/2: L{9) = L„ 



29 (B4) 



Since each bubble edge exerts a pulling force of magni- 
tude A along the outward normal vector of the obstacle 
boundary, the network contribution to the drag is 

To compute this integral, we introduce two dimen- 
sionless variables, both functions of A and R: 



R 

I' 



2R 



(B6) 
(B7) 



The physical meaning of s is equivalent to the liquid 
fraction, since 



$ = (2\/3-7r)£ 



(B8) 



On the other hand, (3 quantifies the amount of 
up/downstream asymmetry, that is, the squashing and 
stretching of bubbles. It increases when <I> (or equiva- 
lently e) decreases (eq. 



). (B9) 



When (f> goes to zero, e goes to zero too, and /3 diverges. 
Using these variables, eq. (jB5[) yields 



F : 



A dp 

LmsLX 



f3 



= arctan (^y/ (3 — ij — 1 



(BIO) 



At high liquid fraction, the force F vanishes when 
P =1, that is (eq. |B9| when: 



(Bll) 



$ = ^^^= 0,086. 

2 + V3 



At low liquid fraction, we develop eq. (jBlOp to leading 
order in /3 and insert the leading order term of eq. (jB9[) 
to obtain eq. ([TU)) . 
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